
%  Modified from the Thermaid code package
udata.ibcp = []; % Empty initialization of internal pressure BC.

%% GRID PARAMETERS ---------------------------------------------------------------------------------%
udata.len     = [130 100];                         % Physical length of the domain in x and y direction [m]
udata.Nf      = [130 100];  	                     % Number of cells in x and y direction
udata.dx      = udata.len./udata.Nf;                 % Cell length [m]

%% SIMULATION PARAMETER FOR TRANSPORT --------------------------------------------------------------%
udata.timeSim  = 30*86400;                     % Total simulation time 100��[s]
udata.dt       = 2*86400;                     % Time step length [s]
udata.tol      = 1.e-3;                            % Tolerance on pressure-heat loop [-]
udata.maxit    = 100;                              % Maximum number of pressure-heat loops to converge

%% INITIAL CONDITIONS--------------------------------------------------------------------------------%   
udata.T0   = 100*ones(udata.Nf(1),udata.Nf(2));    % Initial matrix temperature [oC]
udata.tmax = 120;                                   % Maximum temperature [oC] for plotting
udata.p0   = zeros(udata.Nf(1),udata.Nf(2))+30*1e6;       % Initial matrix pressure [Pa]

%% BC FLUID ----------------------------------------------------------------------------------------%
udata.ibcs = zeros(2*sum(udata.Nf),1)+0;             % Type 0:Neumann(N); 1:Dirichlet(D)
udata.Fix  = zeros(2*sum(udata.Nf),1);%30*1e6;             % Value N [m2/s] (inflow>0); D [Pa]   

udata.ibcs(round(1:2*udata.Nf(2))) = 1;
udata.Fix(round(1:2*udata.Nf(2))) = 30*1e6;

% udata.ibcs(1:udata.Nf(2)) = 1;
% udata.Fix(1:udata.Nf(2)) = 2.5e7;
% 
% udata.ibcs(udata.Nf(2)+1:2*udata.Nf(2)) = 1;
% udata.Fix(udata.Nf(2)+1:2*udata.Nf(2)) = 0.0;

%% BC TRANSPORT ------------------------------------------------------------------------------------%
udata.flagHeatTransport = 1;                       % If 1 - heat transport is used
udata.FixT     = 100*ones(2*sum(udata.Nf),1);      % Temperature at the boundary [oC]
% udata.FixT(round(udata.Nf(2)/2):udata.Nf(2)) = 80;

%% GRAVITY---------------------------------------------------------------------------------%
udata.gravity   = 0;                               % Gravity acceleration in y-direction [m/s2]

%% PERMEABILITY ------------------------------------------------------------------------------------%
udata.K       = ones(udata.Nf(1),udata.Nf(2))*1e-15;% Permeability field [m2]

%% Porosity ----------------------------------------------------------------------------------------%
udata.phi     = ones(udata.Nf(1),udata.Nf(2))*0.3;% Porosity field [-]

%% Fluid density ------------------------------------------------------------------------%
udata.const_density = 1;                           % If 1 - constant fluid density is used
if(udata.const_density)
    udata.density_l  = 1000*ones(udata.Nf);        % Density of the fluid [kg/m3]
end

%% Fluid viscosity ------------------------------------------------------------------------%
udata.const_viscosity = 1;                         % If 1 - constant fluid viscosity is used
if(udata.const_viscosity)
    udata.viscosity   = 1e-3*ones(udata.Nf);       % Viscosity of the fluid [Pa*s]
end

%% ROCK DENSITY ------------------------------------------------------------------------%
udata.density_s  = 2700*ones(udata.Nf);            % Density of the rock [kg/m3]

%% MECHANIC PROPERTIES
udata.flagIncompressible = 1;                      % If 1 - incompressible fluids are used
udata.compressibility_l = 1e-59;                 % Compressibility of the fluid [1/Pa]
udata.compressibility_s = 1e-59;                 % Compressibility of the rock  [1/Pa]
udata.shear_modulus = 29e9;                        % Shear modulus of the rock [Pa]
udata.poisson_ratio = 0.25;                        % Poisson ratio of the rock [-]
udata.K_enh = 1e2;				                   % Permeability enhancement factor [-]
udata.therm_exp_coeff = 7.9e-6;                    % Thermal expansion coefficient of the rock matrix [-]
udata.flagFracStability = 0;                       % If 0 the remaining parameters in this section are not used
                                                   % TR - Trend  (from N) / PL - Plunge (from horizontal)
udata.stress_trend  = [0 90];                      % [TR(S1) TR(S3)] [°]
udata.stress_plunge = 90;                          % PL(S1)          [°]
udata.use_thermal_stress = 0;                      % Boolean for thermal stress in calculations.  If 0 the contribution of thermal stress is not calculated.

%% THERMAL DIFFUSION ------------------------------------------------------------------------%
udata.lambda_l = 0.5;                              % Thermal conductivity of the fluid [W/(m*K)]
udata.lambda_s = 2.0;                              % Thermal conductivity of the rock [W/(m*K)]
udata.cp_l = 4000;                                 % Specific heat capacity of the fluid [J/(kg*K)]
udata.cp_s = 1000;                                 % Specific heat capacity of the rock [J/(kg*K)]
udata.ibcD    = zeros(2*sum(udata.Nf),1);          % 1 -> Diffusion on boundary cells
%udata.ibcD(round(udata.Nf(2)/2):udata.Nf(2)) = 0;
